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Abstract 

Common time-explicit numerical methods for kinetic simulations of plasmas in the low- 
collisions limit fall into two classes of algorithms: momentum conserving [also known as 
Particle-In-Cell (PIC)] and energy conserving. Each has certain drawbacks. The PIC 
algorithm does not conserve total energy, which may lead to spurious numerical heating 
(grid heating). Its overall accuracy is at most second due to the nature of the force inter- 
polation between grid and particle position. Energy-conserving algorithms do not exhibit 
grid heating, but because their formulation uses potentials, computationally undesirable 
matrix inversions may be necessary. In addition, compared to PIC algorithms for the 
same accuracy, these algorithms have higher numerical noise due to the restricted choice 
of particle shapes. Here we formulate time-explicit, finite-size particle algorithms using 
particular reductions of the particle distribution function. These reductions are used in 
two variational principles, a Lagrangian-based and a Hamiltonian-based in conjunction 
with a non-canonical Poisson bracket. The Lagrangian formulations here generalize pre- 
vious such formulations. The Hamiltonian formulation is presented here for the first time. 
Many drawbacks of the two classes of particle methods are mitigated. For example, re- 
strictions on particle shapes are relaxed in energy conserving algorithms, which allows to 
decrease the numerical noise in these methods. The Hamiltonian formulation of particle 
algorithms is done in terms of fields instead of potentials, thus avoiding solving Poisson's 
equation. An algorithm that conserves both energy and momentum is presented. Other 
features of the algorithms include a natural way to perform coordinate transformations, 
the use of various time integrating methods, and the ability to increase the overall ac- 
curacy beyond second order, including all generalizations. For simplicity, we restrict our 
discussion to one-dimensional, non-relativistic, unmagnetized, electrostatic plasmas. 
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1. Introduction 

Particle-based simulation methods have a long and successful history in plasma physics. 
The original proposal by Hockney [T] , based on a method developed for fluid simulations 
[2], used computational macro-particles moving in a self-consistent, mean field. Fields 
were approximated on a spatial grid and interpolated to the particle position to determine 
the Lorentz force. The essential physics was successfully captured but the simulations 
suffered from high numerical noise as (^-functions were used to represent the macro- 
particles. It was later realized [3], that 5-function particles used in this way can lead 
to numerical instability [3]. A significant improvement was achieved by allowing the 
macro-particles to have a finite spatial extent [Jj ; these improved schemes originated the 
class of methods now known as Particle-In-Cell (PIC) algorithms. The PIC algorithm 
was first used to model plasmas with negligible collisions but later Coulomb collisions 
and atomic physics were included with techniques based on the Monte Carlo method [6] . 
PIC simulations are now widely used to study far ranging plasma systems including 
modeling laser-plasma interactions [THllj. Z-pinches fT^, astrophysical and magnetized 
plasmas |13H15] , plasma discharges and low temperature plasma processing [TH] , and 
numerous other applications. 

PIC methods have subsequently undergone a significant development. Research on 
improvements in reliability and stability lead to the discovery of non-physical, purely 
numerical artifacts. The PIC algorithm does not conserve total energy exactly (even in 
the absence temporal discretization), which leads to a surprising phenomenon known as 
"grid heating" [T71 [TH]. (See Ref. [15] for a overview of grid heating.) Grid heating is 
attributed to a kinetic instability where sub-grid structures are aliased to low frequency 
modes (due to finite grid resolution) . Typically this instability saturates once the plasma 
has heated to the point where the Debye length is on the order of the grid spacing [T7]. 
More recently, it has been shown that the choice of particle shape (the spline used for 
current and charge deposition and force interpolation) can lead to unphysical effects, 
especially when considering threshold phenomena such as self-trapping in a laser-plasma 
accelerator jl9] . A further limitation of the PIC algorithm is that its overall accuracy is 
at most second order in both space and time. This is due to the interpolation (typically 
splines) of quantities between the continuous particle position and the spatial grid. 

In an attempt to correct for these deficiencies of the standard PIC algorithm, energy- 
conserving particle algorithms were devised |20j . While strict energy conservation elimi- 
nated the grid heating instability, these algorithms had their own drawbacks. They did 
not seem to have the same flexibility with respect to a choice of particle shapes as PIC al- 
gorithms, which in turn affected the level of numerical noise; i.e., for the same numerical 
accuracy, PIC algorithms had lower noise levels. Energy conserving algorithms also may 
require mass-matrix inversions, which is avoided by the field-based formulations of PIC. 
As a result, energy-conserving algorithms did not become as popular as PIC algorithms. 

One major difference between PIC and energy-conserving algorithms lies in the way 
each is formulated. In PIC algorithms [3T1 [52] , relations between the discretized electric 
(vector) potential, electric (magnetic) field, charge deposition, and current deposition 
were obtained by discretizing the corresponding continuous relations and equations. In 
this process, critical terms of the order of the accuracy of discretization are dropped (as 
is justified in an asymptotic procedure), but which led to the loss of energy conservation 
(and possibly violation of other conservation laws). In comparison, energy conserving 
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algorithms were derived from a variational principle [lUJ [53], using the fact that the 
Vlasov equation could be obtained from Low's Lagrangian [21] . A number of benefits in 
using variational principle was pointed out: basic properties of the original system were 
retained in the reduced system; a natural way of making coordinate transformations was 
provided; and use of high accuracy space and time solvers was possible. 

The goal of this paper is to generalize previous variational formulations and to offer 
a new formulation of energy conserving algorithms based on the Hamiltonian and the 
non-canonical Poisson bracket proposed by Morrison [55H27j . Our approach is to use 
particular reductions of the distribution function to a finite collection of terms as well 
as particular reductions of the continuous fields to finite number of degrees-of-freedom, 
either in the Lagrangian or in the Hamiltonian and Poisson bracket. As a result of our 
general method, we show how to avoid many of the previous drawbacks and deficiencies 
of both PIC and the energy-conserving methods. 

In addition to the energy-conserving property of all algorithms in this work, we: 
(i) show that particle shapes in energy conserving algorithms can be chosen with more 
freedom instead of being a delta- function in space. In fact, the shape of an extended space 
particle has very few physical constraints; the shape may be symmetric about its centroid 
or may exploit the spatial symmetry of a particular physical problem, e.g., systems with 
azimuthal symmetry or symmetries in the gyrokinetic approximation |281 129] ; (ii) relax 
the method of time integration of the equations of motion (most often leapfrog previously) 
allowing for higher than second order accuracy. The choice of a time integrator becomes 
limited only by numerical stability. In this paper we emphasize the spatial discretization 
and leave time continuous, which allows for convenient formulations of time-explicit 
schemes. We prove conservation laws with continuous time and only in the last step do we 
choose a particular (explicit) time advancing scheme; (iii) show how to reduce continuous 
quantities with either grid-based reduction (using finite differences) or truncated bases; 
previous authors have only used truncated bases, which for the case of finite elements may 
necessitate mass matrix inversions. By using grid-based reduction, mass matrices do not 
appear (for an example, see Appendix B ), which may be computationally advantageous; 
(iv) derive formulations in terms of fields that are based on the Hamiltonian and a non- 
canonical Poisson bracket. Previously only potential-based energy-conserving algorithms 
have been derived from a variational principle [201 ES] . Using field-based formulations 
eliminates the need for solving Poisson's equation; and (v) derive a particle algorithm 
that conserves both total energy and total momentum. The arguments leading to this 
algorithm demonstrate the usefulness of the variational approach and exploit the relations 
between conservation laws and symmetries of the Lagrangian. Previous such models were 
derived by assuming a delta function for the particle shape [301 131j . Throughout, we 
restrict our discussion to the case of a one-dimensional, nonrelativistic, unmagnetized, 
electrostatic plasma. This is done purely to streamline the discussion, elucidating the 
central ideas. The extension of these concepts to the fully relativistic, electromagnetic 
case, while involving numerous technical details, is largely straightforward and will be 
the subject of a future publication. 

It has been emphasized before [20l [23] that variational formulations naturally lead 
to higher overall accuracy particle algorithms (i.e., both in space and time). Here we 
show that this remains true with all of the above relaxed conditions. We show that 
force interpolation, field integrators, and time integrators may be chosen to increase the 
accuracy of a particle method beyond second order. In the course of all derivations, we 
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point out where a certain property is being relaxed or is being lost. 

Time-implicit formulations of particle algorithms have significant attraction in sim- 
ulating problems where long time evolution is necessary. Recently, authors have been 
successfully reformulated the PIC algorithm in terms of implicit time integration with 
the added benefit of energy conservation [32[ I33j . However, these formulations do not 
offer a general derivation and so is unclear how they can be extended. For example, 
these formulations use the Crank-Nicholson time integrator and charge-conserving parti- 
cle shapes but it is an open question how to extend these methods to use more accurate 
time integrators. 

While continuous-time formulations are the focus of this paper, we note that one may 
discretize the action principle in both space and time. In the context of particle methods 
for plasma simulations, this was first considered by Eastwood [23 . While Eastwood's 
method conserved energy exactly, it did so at the expense of an implicit time-advance. 
When the time-advance was altered to be fully explicit, exact energy conservation was 
lost. The notion of performing the temporal discretization in the action has been studied 
extensively by Marsden and co-workers (see for example Refs [34] and [35]). There are 
a number of attractive features of this approach and it is a natural extension of the 
methods presented herein; this will be the subject of future work. 

The paper is organized as follows. Section [2] is devoted to deriving algorithms based 
on a Lagrangian formulation of the Vlasov-Poisson system. This section presents the 
finite differencing formulation. In section 2.2 the particle models are derived from a La- 
grangian formulation in terms of truncated bases. The important model which conserves 
both momentum and energy is derived there. In these derivations the fields are described 
in terms of electrostatic potential. Section [3] presents the Hamiltonian derivation of par- 
ticle models using truncated basis and a reduction of the non-canonical Poisson bracket. 
The equations of motion are formulated in terms of electric field. Section [4] illustrates 
properties of the derived particle models with numerical examples. Conclusions are in 
section [5] [Appendix A| gives many examples of charge deposition rules, while [5"ppcndix| 
|B] presents a hybrid cold fiuid-kinetic particle model from a Lagrangian starting point. It 
demonstrates our general method with a different reduction of the particle distribution 
function. It also illustrates how the mass matrix and its inverse may be avoided by the 
use of grid-based reduction of the continuous quantities. 



2. Lagrangian formulation 

A plasma with negligible collisions is well described by a single-particle phase-space 
distribution function, /, whose phase-space evolution is governed by the Vlasov equa- 
tion [3^ 

— +u— + — i;— =0 (1) 
dt dx nis dv ' 

where E = —Vip is the electric field, ip is the electric potential, and Qs are the species 
mass and charge. For an initial phase-space distribution fo(x,v), the distribution at any 
later time is given by 

f{x,v,t) = foii.i), (2) 

where x{t;x,v) and v{t;x,v) = dx{t;x,v)/dt are the macro-particle trajectories with 
initial conditions x and v: x{0;x,v) = x and v{0;x,v) = v. The particle trajectories 
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correspond to characteristics of the Vlasov equation and ^ is simply the statement 
that the distribution function is constant on the characteristics. Vlasov dynamics can be 
obtained from the Lagrangian [231 ISZl 



£ = ^ jdxdv fo{x,v) 



dx{t] x, ■0)"' ^ 
dt 



-qs jdxdv foix,v)ip{x{t;x,v),t) + ^ Jdx [V(p{x)f , (3) 

where x{t; x, v) and f{x) are to be varied independently. (Here we consider a single- 
species plasma but the extension to multiple species is obvious.) In the usual way, 
demanding the action be stationary with respect to variations of the dynamical variables 
leads to the equations of motion. Variation with respect to particle positions gives 

rUsX ^ -qsV^ , (4) 

while variation with respect to the potential gives 

V^</3 = -47rgs jdv f{x,v,t) = -At:p{x). (5) 



We have used [37] 

dx dvfo{x,v) — dx dvf{x,v,t) (6) 

in ([5]) and we have assumed either periodic boundary conditions or an infinite system 
to allow surface terms to be dropped. Note that (|6]) is a statement of particle number 
conservation and is equivalent to Gardner's restacking theorem |39j . 

The basic idea of Lagrangian macro-particle methods lies in the representation of the 
full distribution function f{x,v,t) as a sum of moving spatial volumes, fa{x,v,t), called 
macro-particles: 

f{x,V,t) = ^fa{x,V,t) 

a 

= ^^«a^[x-ea(iPb-€a(t)]. (7) 

a 

The choice of a delta function in velocity space is not essential but avoids the necessity 
to track stretching phase space volumes, which is why we adhere to it. In Eq. ([?]), Wa 
are constant weights and the function S is the fixed spatial extent of the computational 
particle (hereafter we use the terms particle, computational particle, and macro-particle 
interchangeably unless otherwise specified) and is normalized as 

dxS[x-^ait)]^l. (8) 

An additional simplification is made by assuming that all particles have the same shape. 
We note that the representation ([T]) is general and independent of whether both electric 
and magnetic fields are present in the system, i.e., it is valid for the general Vlasov- 
Maxwell system. For clarity of the presentation, in this paper we consider only electro- 
static, non-relativistic models. Electromagnetic and relativistic models can be derived 
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similarly to this presentation and will be presented in a future publication. We view ([T]) 
as a particular reduction of the particle distribution function. [Appendix B| gives another 
example of such a reduction. 

Substituting our form of the distribution function, into the Lagrangian and again 
using (|6|, we obtain a reduced Lagrangian 

m • /■ If 

^^Wail- qs^Wa jdxS{x-ia) ip[x) + — dx Ifif 

a—l Q— 1 ^ ^ 

— -^kin ~^ -^int ~t~ 'Afield i 

where 

Uin^'^Y.W^il (10) 

a=l 

Cint= Jdx S{X - (,a)f{x) , (11) 
Acid = ^ Jdx {Vip f. (12) 

Although we have replaced a continuum of particles with labels x and v by Np macro- 
particles, we still have an infinite degree-of-freedom system due to the presence of the con- 
tinuous field (fi. The equations of motion are obtained from (|9| by considering variations 
of the particle position and of the potential. For the particles, the usual Euler-Lagrange 
equation 

±dC_dC 



L = -— [dx ^ip{x) ^ -— ldxS[x^^a.{t)]Vip. (14) 



gives 

nis J d^a ^ ^ ^ m 
Since the potential is a field, the Euler-Lagrange equation for the potential is 

^ = 0, (15) 
d(p 

where 5 /Sip denotes a functional derivative. Then 

VV = -4vrg,^u;„5[a;-^„(t)]. (16) 



Note that the factor qs/rus appearing in (14) is the physical charge to mass ratio of 



the plasma species. It is not necessary to make the ad-hoc assumption that the macro- 
particle have the same charge to mass ratio as the plasma species, this is a consequence 
of the phase-space decomposition ([7]). Furthermore, the second form of the force in (14) 
may clearly be interpreted as the electric field averaged over the particle shape. 

The substitution of ([t]) into ([s]) is equivalent to a choice of a trial function for / 
that depends on a number of parameters, which in our case are the particle positions 
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and velocities. The values of these parameters are obtained by solving the equations 
resulting from the variation (13). Other choices of trial functions may lead to models 
without particles at all [JHl EI] • 

A significant advantage of the variational formulation is the connection between sym- 
metries and conservation laws as embodied in Noether's theorem [42 . Our introduction 
of macro-particles through ([t]) neither results in explicit time-dependence in the La- 
grangian nor breaks translational invariancc of the Lagrangian, thus we should expect 
the equations of motion ( |14[ ) and ( 16 ) to exactly conserve both energy and momentum. 
The total energy of the system is the sum of macro-particle kinetic energy and field 
energy, 

W^=^E^"^' + ^ dx{Vv)\ (17) 

Q — 1 

Using the equations of motion, it is straightforward to see that W is an invariant: 

M 

dW ^ ■ •• If, 9 ^2 

— = rris 2_^Wa^aia - ^ dx (f — V if 

a— 1 

g iVp 

= -qs ^ Wa ia dx S{X - Ca)V(^ + Ux V WaS{x - ^a) 

a=l a=l 

= ~qs '^Waia / dx S{X - ^a)V(y9 - Qs'^Wata dx (f -^S{x - 
Q=l a=l ■' ^ 

r . /■ 

= -qs '^Waia / dx S{X - ^a)V(y9 + qs'^Waia dx S{x- £,a) V(p 
a=l a=l 

= 0. (18) 
The total momentum of the system is simply 



P = ms'^Waia, (19) 

a=l 

since the electrostatic field carries no momentum (the Poynting vector is zero in the 
electrostatic approximation). Now 



dP 
~dt 



-^s^Wq jdx S{x - ia)^(p 

dx V'^ipVip 
dxV {Vip f 



1 



47r 
1 

8^ 

0, (20) 



where we have used (16 1. 
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At this point in our reduction, we have a finite number of macro particles represent- 
ing the plasma but a continuous field for the potential. Here one must provide some 



approximate or exact solution of (16), which is then used to integrate the macro-particle 
equations of motion. One possibility would be to use methods based on evaluating the 
Green's function, constructing ip as the superposition of the potentials due to each macro- 
particle. Even though the macro-particles interact via the mean field, computation of Lp 
using Green's functions scales as OiN"^) and is thus limited to relatively small systems. A 
more computationally advantageous alternative is to introduce a discrete representation 
for the potential. There are two general approaches: using a spatial grid, approximating 
the potential by its value at the grid point, or using a truncated set of (local or global) 
basis functions and representing the potential by its projection onto the basis. 



The interaction term in the Lagrangian, (11), provides both the force in (14) as well 



as the charge density in (16). Of course, this will continue to be the case when the 
continuous potential is replaced by a discrete approximation. It will be necessary to 
approximate >Cint consistently with the choice of the discrete potential but this single 
approximation ultimately yields both the force term in the equation of motion and 
the charge density in a discrete analogue of Poisson's equation. Thus we are guaranteed 
that these terms are consistently approximated. 

2.1. Discretization using a spatial grid 

We assume a fixed spatial grid Xi with i e [1, A"g] and grid spacing h with (pi being 
the numerical approximation to (p{xi). We must now approximate two terms in ([9|, £int 
and £fioid- The interaction term requires knowledge of tp between the grid-points; so 
some manner of interpolation is required. Finite elements |43j offer a consistent way to 
perform such interpolations to any accuracy. Let ^^(a;), i = 1, . . . ,Ng be finite-element 
basis of some order. We interpolate p between the grid points by 

(^(x) = ^^,*,(a;) (21) 



and thus (11) becomes 

/- ■ y /• - ■ y 

dx S{X - £_a)'P{x) ^'^Pi dx S{X - £,a)'^i{x) ^'^PiPi{^a) , (22) 



1=1 



where 

= JdxSix-^^)^,{x) (23) 

is the effective (projected) shape of the particle. Note the expression for pi can be 
computed analytically since the function S is known. If ^',(2;) are constructed from 
Lagrange polynomials, then ^^ii^) = 1 a-nd 

= ldxSix-^^) = l. (24) 

i—l i—1 

This property means that the total charge deposited on the grid and at any instant of 
time is constant. 



It remains to approximate >Cfioid in terms of (fi . This can be approached in two ways 
which give roughly equivalent results. We can use (21) to write the integral in (12 1 as 



Defining 



we have 



hKj^ = I dx 



d'^i{x) d^j{x) 



dx 



dx 



jdx {S/(p f ^ "h'Y^ ip^tfjKij . 



(25) 



(26) 



(27) 



Alternatively, after integrating by parts in Afield, we can approximate the integral as 



dx (V93) 



J dx tp{x)V'^ ip{x) « ~h 



dx"^ 



(28) 



While this appears to simply be using the trapezoidal rule to evaluate the integrand, 
with either periodic boundary conditions or an infinite domain, this approximation has 
spectral accuracy, that is, all modes supported by the grid are integrated exactly. We 
complete the approximation by choosing a finite-difference representation for the second 
derivative. Regardless of details of the finite-difference approximation, it can always be 
expressed as 

d^if 



dx^ 

for some integer a. Thus we have 



J = l 



dx (V^pY 



-h ^ ipiK^-jipj , 



(29) 



(30) 



which has the same form as (27). Notice that while Kij is always symmetric [cf. (26)], 
this need not be true for Ki^. 



We now arrive at the finite degree-of-frecdom Lagrangian 
m h 



(31) 



where we take either )Cij — Kij or /C^ — Kij . The dynamical equations are obtained by 
demanding that the action be stationary with respect to variations in both S.a and ipi. 
Taking these variations yields 



N„ 



Co 



(32) 



and 

'y'j = -47r ^ ^ Pi ) . (33) 

3 = 1 ^ a=l 



The reason we choose to integrate by parts in ([30|) is now clear: by dong so, we are 
able to directly specify the difference method for the second derivative that appears in 
Poisson's equation, (33). 



Discretizing (p{x) in (17) in the same manner as in the Lagrangian, we have 

WL = '^j:w^il-^Y.^^^^^^r (34) 



Using the equations of motion, wc find 

oWl ; 1- h y-^ ^ dip. 



N„ N„ „ N„ N, 



-9; 



dpi , dpi{ia 



a—l i—1 i—1 a—l 

\ ^ \ ^ Cpi ■ \ ^ \ ^ Opi ■ 

= -9s 2^ 2^ ia+ q8 2_^2^Wa 

= . (35) 

Introducing a spatial grid does not affect energy conservation. This is expected since 
the spatial discretization of ip does not introduce explicit time-dependence into the La- 
grangian. An immediate advantage of the variational approach is that models derived 
in this way are automatically free of grid heating. A side effect of introducing a spatial 
grid is that it breaks the translation invariance of C and consequently total momentum 
is no longer exactly conserved; see section 2.2 for a more complete discussion. 

We conclude this section by providing a concrete example of this procedure to derive 
a model that is second-order accurate in h. For simplicity, we consider the case of 
a charge-neutral electron plasma with an immobile ionic background and a spatially 
periodic domain. Since this system has no ion dynamics, we can forgo summing over 
species and simply introduce the ion density into the Lagrangian 

^ ^^Wo^il- qe^Yl '^o.PtiCa) 'Pt'Yl '^^ ~ ^ X! ftK-'ijfj, (36) 

a — l a—l i—1 i—1 ^J — 1 

where 

pT'' ^qjdxn''°^\x)^,{x), (37) 

with ri''™' (x) being the given ion density. Linear finite elements yield second order accu- 
racy interpolation [321 (see Figure [T]): 

{\x — Xi\ 
1 ^ x,_i<x<x,+i, ^3g^ 
otherwise. 
10 



To determine pi we need to specify S{x). Regardless of the choice of S, the accuracy of 
the interpolation will be second order due to our basis choice. The choice of S affects 
the quality of our approximation through the extent to which ([t]) is a good ansatz but 
has no influence on the formal order of the model. Arguably the simplest choice for S is 
a top-hat cell-wide function: 



otherwise. 



(39) 



While we have chosen S to be exactly one grid-cell wide, this is by no means essential. 
The choice of support of S is completely independent of the grid spacing; the particular 



choice in (39), allows us to make connection with the usual PIC particle shapes (see 
Figure |A.9|and Table IATI) . 




Figure 1: Linear finite element basis functions. The basis functions are identified witfi tfie grid-point 
at wfiicfi it takes on tfie value 1; e.g., '^i is the tent function with support [x'i_i, Xi_|_i]. 



We now use (23) to determine the grid charge deposition. With this S, for any 



there are only three values of i for which Pi{£,a) 0. Take Xk to be the grid point nearest 



and let A = {(.a — Xk)/h. Clearly |A| < 1/2. It is straightforward to evaluate (23 1 to 
obtain 



Pk-i 
Pk 



2 V 2 



(40) 



This is equivalent to the charge deposition obtained from the usual PIC quadratic particle 
shape [21]. It is possible to recover all of the usual smooth particle shapes. For example, 
taking 



1 







otherwise, 



11 



we obtain 



Pk-2 = 



Pk-i 



Pk+i 



24 

19 

96 



A- 



11 
24 



A- 



115 5 ^2 
=192^8^ 



19 
96 



11 
24 



A^ 



1 



- A3 - ^ A* 
6 6 



6 



6 



;a^ 



(42) 



Pk+2 - 



which is equivalent to the usual quartic charge deposition rule. We take up the matter 
of particle shapes in some detail in [Appendix A[ In particular, we demonstrate a cubic 
Pk spanning only three grid-points. 



All that remains is to approximate ( 12 ) to second-order accuracy. Evaluating ( 26 ) 
'^'^^ is straightforward. For any i, we can see that '^f^'{x) has a 
'^''(x) only for j — i±l and thus we have 



for our linear basis 
non-zero overlap with 4'^ 




(43) 



i±l. 



Alternatively, we can use finite difference approximations and evaluate ( 12 1 using 



(30). Since we are considering a periodic domain, it is reasonable to choose a central 



difference approximation 



0{h^), 



which gives 



(44) 



(45) 



Linear finite-elements give the same approximation to ( 12 ) as taking second-order central 



differences. This is essentially a coincidence; higher-order finite element bases (quadratic, 
cubic, etc.) do not yield expressions for Kij that can be equated to conventional differ- 
encing schemes. Nothing, other than the symmetry of the problem, forces us to choose 
central differences; any second-order approximation would suffice. Note, using an ex- 
pression for Kij that is accurate beyond second-order will not increase the overall spatial 
order of the method unless a corresponding more accurate interpolation scheme is used 



to evaluate ( 23 ) 



The macro-particle equation of motion, which is not alerted by the specifics of the 



spatial discretization, remains as in (32). For a uniform ion background, rig™ , we have 



(Ion) 



12 



(Ion) 



h. 



(46) 



For completeness we restate Poisson's equation including the ionic background and our 
approximation for ICij 

- 2(pi + ip^+i) = -Att -^^Wa Ptiia) - 47rgi no™' . (47) 

2.2. Discretization using a truncated basis and the question of momentum conservation 

In the previous section the interaction and field parts of the Lagrangian (|9| were 
reduced from an infinite to a finite degree of freedom quantities by discretizing on a grid. 
We noted that the introduction of a spatial grid breaks the translational invariance of 
£, which leads to loss of momentum conservation in the reduced system. In this section 
we consider a reduction using a truncated global basis and investigate the question of 
momentum conservation in this case. We show that replacing the continuous potential 
by a finite collection of projections onto a truncated basis can result in a discrete system 
that retains translation invariance. (Of course, if the basis is not truncated, which is not 
useful from a computational perspective, then we would expect translation invariance to 
be maintained for any complete basis.) 

Let $m(a;), m = 1,...,M be the first M elements of an orthonornial basis. We 
approximate the potential as 

M 

where 



iPm= jdx^i,{x)if(x) (49) 
and '^m{x) is the dual to $,„(x) satisfying 

\dx^,n{x)^i{x) ^ 5ran, (50) 



for m,n = 1, . . . ,M. The accuracy of this approximation will depend on the number 
of elements kept and the convergence properties of the basis. For a truncated system 
to possess translation invariance requires that the basis have certain properties. These 
properties are made evident by shifting the origin by an amount 5x while leaving the 
physical system unchanged and requiring this transformation to be a symmetry of C '42j . 
The kinetic term, /Ikin, is obviously translation invariant under such a shift since 5x is 
time independent. Consider £int. Let x' be the new coordinate, with x — x' + Sx. The 
particle coordinates and potential relative to x' are 

— ~ (51) 

ifiix') — ip{x' + Sx) . 

The symmetry condition is then 

-Cint Ca] = -Cint |a] • (52) 
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Introducing our basis expansion into £int and suppressing prime symbols, we have 

Np M 



a—1 m — 1 

Np M 

= -Qs X] X ^™ / S{x- Sx)^m{x) 
a=l m—1 

Np M „ 
= -Qs X] X "^^ /^^ "^(^ ~ £.a)^rn{x ~ 5x) 

^ 1 ™ 1 



(53) 



(Pm ^ dx ^ln{x) ip{x) = dx ^ln{x) ip{x + Sx) ^ dx ^l^{x - Sx) ip{x) . (54) 



and 



Combining these expressions and expanding to lowest order in Sx, we have 

Np M 



Np M ( „ 

a=\m=\ \ 



dxS{x- £,a)^r,i{x) 



Sx 



dx f{x) 



d^lix) 



dx 



dx S{X - ^a)^7nix) 



+ fdx ip{x)^l,{x) fdx S{x- £,a) ^^"^^^^ 



Np M 



dx ^{x) 



dx 



dx S{X - ^a)^rn{x) 



dx ipix)<i>lix) fdx Six -U^^^^ 



(55) 



Our symmetry condition requires that the term multiplying Sx in (551 vanish. For the 
symmetry to exist independent of the particle shape 5* and the details of the potential, 
this term must vanish for each m. We are led to the condition 



d^^ (x) 

= ±a{m) ^,n(x) and — ^j^-^ = =Fa("i)<l'j„(a;). 
dx dx 



d^m{x 



(56) 



On a finite domain with periodic boundary conditions this condition is satisfied by the 
discrete Fourier basis; we are not aware of any other discrete basis that fulfills (56) on 
either a finite or infinite domain. Using the discrete Fourier basis, it is straightforward 
to show that £fioid is also translation invariant. 

We now specialize our discussion to the case of a truncated Fourier basis. Let 



I 



ikx 



(57) 



L 
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where k = 2rmr/L, m = 0,±1, . . . ,±M and L is the domain size. With this basis, the 
interaction term becomes 



iO-int = -qs^^Wa(pk / dx S{X - ^a)^ki3 



a—1 k 

-qsL^^Watpk P*k{£.a) 
a—1 k 







(58) 



where 



and 



Pk{^a)= / dxS{x-^a)<^l{x), 



ipk= dxip{x)<^l{x), 



(59) 
(60) 



and we have used the relation $fc(x)/L ~ [^^.{x)]* . We also need to evaluate (12): 



d^kjx) d<S>k'{x) 
dx dx 



= -—S^ kk'ifikVk' / dx^kix)^k'{x) 

VfcfcVfc'^fe' / dx<^>k{x)['i>l,{x)]* 
I. Jo 



L 



(61) 



L 

'8^ 



^kk'ifikVk' f dx<^k{x)<^i 
k,k' Jo 



where, since is real, ip-k = Va;- Finally we arrive at the discrete form of the Lagrangian 



a—1 A; 



47r ■ 



(62) 



fe>0 



To obtain the equations of motion, we require the action to be stationary with respect 
to variations of and Lpk (since (pk and are not independent, we need only consider 
variations of ^Pk)- The equation of motion are 



E 

k 
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fk 



(63) 



and 



Using ( 59 1 and ( 57 ) it is easy to show 
allowing us to write the equation of motion as 



qsL 



qsL 



k>() 



^2/cIni[p*(C„)^fc] 



fc>o 

The spatial charge density associated with a single particle is 

qs^Pk{^a)'^k{x) 
k 

and the corresponding total charge is 



u Jo ,. 



qsLpQ{^a) 



qsL f dx S{x - Ca) Y 
Jo ^ 



(64) 



(65) 



(66) 



(67) 



(68) 



Thus, regardless of the number of modes retained, the charge associated with each particle 
remains qg- 

The energy of this system is 



Using the equations of motion we have 

^ 7/)„. ^„ -I- 

47r 



(69) 



fe>0 



dt 



k>0 



iqsL^^w^kia {p*k Vk - Pk Vl) - iqsL^^kia {PkVl " Pk Vk) 
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a=l fe>0 



fc>0 Q = l 



(70) 



where we have used (64 1 and (65) to find ipk. From (19) we have 



~dt 



Q=l fe>0 

L 



An 



fc^ ifi - fk ^i) 



k>0 



0. 



(71) 



Thus the model using a truncated Fourier basis conserves both energy and momentum. 
This is as expected since the spatial discretization does not introduce time-dependence 
into the Lagrangian and the basis was specifically constructed to maintain spatial trans- 
lation invariance. 

Consider the same system as in Section 2.1 with S given by (39) where h is an 
independent parameter. Now ( 59 1 becomes 

1 



PkiCa) 



L 



e sine 



ilkh), 



(72) 



where sincx = sin(a:;)/a::. If, as before, we take a quasi-neutral plasma with a uniform 
ion density, then the ions only contribute to Poisson's equation for k = 0. Further, we 
see that fc = does not contribute to ^q, and thus we are free to take (po = 0. With this 
form of pk , the potential becomes 



k^L 



sine (ifc/i) ^ Wq. e fc>0 



(73) 



and the particle equation of motion becomes 

Is ; ..-..^^ r 1 r jk^ 



Co 



^(-i fc) sine (ifcft.) [e*''^° ipk 

k>0 

J2 sine (ifc/i) [e'''^-Ek + e-'^^" El] 

fc>0 

^sinc(ifc/i)i;fce''=«° 



(74) 



where Ek = —iktpk- For the complete basis, by the convolution theorem, ( 74 ) is identical 
to (14). While, due to the truncation, the convolution theorem does not apply, we may 
still interpret the force in (74) as sampling the electric field over the effective spatial 
extent of the particle. 

In this subsection we have derived a particle algorithm that preserves the time and 
space translational invariance of the Lagrangian and thus conserves both energy and 
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momentum exactly. Since the use of grid in the reduction violates the spatial translational 
invariance, we were lead to use a continuous basis. In the course of the derivation, we 
found that only one such basis exists, a (possibly truncated) Fourier basis. One may argue 
that the use of a Fourier limits the applicability of this algorithm; for example, restricting 
to systems with periodic boundary conditions or being unsuitable for large-scale parallel 
simulations. This result establishes, however, that simultaneous conservation of energy 
and momentum is indeed possible. 



3. Noncanonical Hamiltonian formulation 



It is well known pSVET] that the Vlasov-Maxwell system possesses a Hamiltonian 
structure in terms of non-canonical field variables. Specializing to the 1-D electrostatic 
case and treating the electric field £^ as a dynamical variable, the Vlasov-Maxwell bracket 
I25l [271 becomes 



{F,G} = Jdxdpf 



6F^ SG 
If' If 



df fSF SG 6G SF 



dp \SE Sf SE Sf 



(75) 



where F and G are any functionals of / and E and [a, b] denotes the usual phase-space 
Poisson bracket: 

\ab]^——-—— (76) 
dx dp dx dp 

The Vlasov equation and the equations for the fields are obtained from this bracket and 
the Hamiltonian 



H = 



1 



2tos 



dxp^ f 



1 

8^ 



dxE^ 



as 



9/ 
dt 

dE 
~dt 



p_ df_ 
m, dx 



dl 
dp 



qs^E, 



= {E, H) = -47r fdp — f^ -47rj. 

J f^s 



(77) 

(78) 
(79) 



Poisson's equation is considered as an initial condition and is satisfied for all time as a 
consequence of ( 79 1 . 

We use a reduction of the distribution function, which is identical to ([t]) but written 
in terms of momentum: 



f{x,p,t) = ^fa{x,p,t) 

a 

= ^ S[X - £,ait)] 6\p - TTait)]. 



(80) 



Consider a single fa- The quantities Wa, £,a, and tTq, which denote the macro-particle 
weight, centroid, and momentum, may be expressed as: 



Wo 



dx dp fa , 
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(81) 



1 

1 



dx dp X fa , 
dx dppfa- 



(82) 
(83) 



Therefore, they may be thought of as functionals of /„. To an^arbitrary functional 
F[f] there exists a corresponding function F{wa,£,a,T^a) such that F(wa, ^a, tTq) = F[f]. 
(Both F and F can also be functionals of E; for the moment, we are only interested in 
their dependence on /). Then a functional derivative of F[f] may be found using the 
chain rule as ^ ^ ^ 

6F _ Swa dF ^ S^a OF ^ Stt^ dF 

Sfa Sfa Wa Sfa S,a 5fa TTq 

Evaluating the functional derivatives of Wq, ^q, and tTq, 



(84) 



SWa 
Sfa 

'5U 



1, 

X ^a 
Wa 

Wa 



we find 



Consider 



SF 

Sfa 



dF 

dWa 



x-^a dF 

Wa d^a 



P ~ TTa dF 
Wa diTa 



'SF 


SG' 




' dF 


SU 


su. 




dWa 



X — £,a dF ^p — Tla dF dG ^ X — ^a dG 
Wa d^a Wa dWa ' dWg Wa dS,a 



dF dG 

d^a d-Ka 

dF dG 

dia dlTa 



X -^a P-T^a 



Wa 

X p 
Wa ' Wa 



dF dG 



dlTa d^a 



p-TTa X-^a 



dF dG 

dlTa d^a 



P X 
W„ ' Wa 



1 



Wf, 



w^ 



dF dG dF dG 

d^a dlTa d-Ka d^a 



F,G 



577 



where 



da db 



da db 



d^a dlTa d^a dlTa ' 



The first terms in ( 75 ) then become 

SF SG 



dx dp fa 



Sfa ' Sfa 



= Jdxdpfa- 
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w^ 



F,G 



577 



(85) 



(86) 



P — TTa dG 
Wa dlTa 



(87) 



Now consider 



I dx dp 



dfo, 5F 5G 



dx dp fa 
dx dp fa 

dx dp fa 



6F d SG 
SEdpJu 

5F d ( dG x-^a dG p-Ha dG 



SE dp \ dwa 
SF 1 dG 



Wa d£,a Wa dlTa 



dG J_ 

dlTa Wa 



6E Wa dlTa 
SF 



dx 



dG 

dlTa 



SF 

ldxS{x-^a)j^, (90) 

where the first fine foUows from integration- by-parts and the fact that SF/SE does not 
have p dependence since is a function of x only. Combining ( 89 ) and ( 90 1 with ( 75 ) 
leads to the bracket: 



{F,G} 



1 



F,G 



4:nqs / dx S{x — 



SG dF SF dG 



SE d-Ka SE d-Ko 



(91) 



We now extend this result to a collection of /„. We treat each fa as a separate species, 
which mandates that the only interaction between the various fa is through the mean 
field. The bracket is thus just the sum of the (91 ) over a: 



{F,G}^Y. 



1 



Wa 



F,G 



-I- Airgs dx S{x - £,a) 



SG dF SF dG 



SE d-Ka SE dlTa 



(92) 



Under our reduction, the Hamiltonian becomes (hereafter we drop the tilde notation 
as it should be clear from the above calculation where a functional derivative or a partial 
derivative is taken) 



a— 1 

and the equations of motion are 

ia — {^a,Et} = — 
Wis 

T^a = {TTa, i?} = j dx S{x - ^a)E{x) 



(93) 

(94) 
(95) 
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dE 

— = {T:a,E} = -47rg3 — S{x - £,a) 

at ^-^ m. 

Q = l 



= -47rgfs ^ 'S'(a; - ^a) = -47r j . 



(96) 



a=l 



Equations (94) and (95 1 are easily seen to be equivalent to (14). Comparing the spatial 
derivative of (96) to the time derivative of (16), we see that (96) and (16) are indeed 
equivalent. As in the Lagrangian case, the reduction from a continuous phase space dis- 
tribution function does not break energy or momentum conservation. In the Hamiltonian 
setting, energy conservation follows from the antisymmetry of the Poisson bracket under 
F ^ G and hence is intrinsic to the theory. 

To complete the reduction to a finite degree-of-freedom model, we represent E using 
a finite, discrete basis, ^Pfe having iVf, elements as 



E{x,t) = Y,E^{t)^r{x). 



i=l 



where 



and 



E,{t)=Y^ dx E{x,t)Mr^H j{a 

3 = 1-^ 



Mij = Jdx 'i'i{x) ^'j(x) . 



(97) 



(98) 



(99) 



When 5'i(a;) are a finite element basis, Mij is called the mass matrix. From (98), we 
have 



6E 



Now, the Ei, through (97), provide a complete characterization of E and thus any func- 
tional of E can be written as a function of the Ei. Consequently 



5 

SE 



SE, d 



SE dE, 



i=i 



d 

dE, 



(101) 



Using this expression, the bracket becomes 

Np Nt Np 

{^'G} = E-[^'G]c.+4^9.EE 

1 "^a ■ -I 1 

a— 1 — l OL—1 



dG dF dF dG 



dEi diTa dEi diTa 



M^^pMc), (102) 



where pj (^q ) is defined by ( 23 ) . The reduction of the bracket is exact in the sense that 
given the representation of / and E, [(80) and (97), respectively] the reduced bracket and 
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full bracket, restricted to functionals of the appropriate form, give the same result. Con- 
sequently, the reduced bracket inherits the Jacobi identity [44] (and all other properties) 
from the full bracket. 

Using ( 97 1 and ( 99 ) , we can write the Hamiltonian as 

a— 1 2J — 1 

The equations of motion are then 

L = — (104) 

TTa = ^ E.i Pi{(,a) (105) 

4=1 

Np Nb Np Nb 

a—1 j—1 ^ a—1 j—1 

(106) 

To make a connection with the model based on finite differences (Sec. [5]), note that 
multiplication by the matrix M is equivalent to performing an integration. For finite 
elements constructed from Lagrange polynomials one may reduce the mass matrix to a 
diagonal form (a procedure known as "lumping") while preserving the accuracy of the 
approximation [45] . If we use linear finite elements on a grid with spacing h, lumping 
the mass matrix gives 

A'U,-^hS,j. (107) 



4. Examples 

In this section we present two examples illustrating some properties of the energy con- 
serving models derived in this paper. We begin with a benchmarking example: the linear 
growth rate of the instability caused by a small electron beam of density ni, propagating 
in a neutralizing background plasma of density uq (beam-plasma instability). For small 
beam-to-plasma density ratio, (rib/no) 1 [more precisely, the parameter (7if,/2no)^/"^ 
must be small in this linear theory] , the linear growth rate of this instability is given by: 

/q / \ 1/3 



All simulations are in dimensionless variables, where the time is measured in units 
of inverse plasma frequency, ujp^, momentum is measured in units of mgC, potential is 
measured in units of nieC^/e, and energy in units of meC^no/ujp (assuming 1-D). In the 
latter notation rrie is the electron mass, e is electron charge. The system is assumed 
to be periodic and its dimensionless size is 27r. In this way, the numerical growth rate 
is dimensionless while the physical growth rate is measured in units of uip. In Fig. [2] 



we show a simulation using the model of Sec. 2.2 (|73|) and (|74|). The simulation was 
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Figure 2: Linear growth and saturati on o f the first four harmonics in the beam-plasma problem com puted 
using the truncated Fourier model, | |73| and | |74| , for nj/no = 10~*. The analytical growth rate, l |l08| , 
gives -yL ~ 0.03190 for k = 1. 



initialized by perturbing the beam density (position of beam particles) at the wavelength 
of the first harmonic, and the velocity of the beam was matched to the plasma wave 
phase velocity, e.g., Wboam = 1- (To initialize the second harmonic, fc = 2, the beam 
velocity would have to be set to Wbeam = 1/2, etc.) The beam-to-plasma density ratio 
for this simulation is 10^^. There are 300 particles for each group of particles, i.e., beam 
electrons, background (plasma) electrons, and plasma ions, as well as 128 Fourier modes. 
The plasma ions neutralize exactly both the beam and the plasma electrons, which is 
achieved by an appropriate choice of particle weight (this assures the potential has zero 
bias). The beam to plasma density ratio was also adjusted by an appropriate choice of 
beam and background particles weight. 

The numerical growth rate of the fundamental harmonic is approximately 0.03132, 
which differs by less than 2% from the theoretical value 7l = 0.03190. (The regions 
where the growth rates are determined are indicated by dots.) Better agreement can be 
achieved for smaller beam-to-plasma density ratios. Also seen from this figure is that 
the next three harmonics grow sequentially as a result of the non-linearity developing in 
the growth of the previous harmonics; i.e., the second harmonic is seeded by the non- 
linearity of the first harmonic when the quadratic term of the field has grown sufficiently, 
etc. In this scenario, linear growth rates of higher harmonics are multiples of the growth 
rate of the first harmonic independent of how well the numerical growth rate agrees with 
formula (108), as long as a clear linear stage exists; this is indeed the case in Fig. [2] 
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Energy and momentum balance are shown in Fig.[3j Momentum is conserved to machine 
precision even in the time-discretized model, while energy conservation depends on the 
time integrator properties and time step At. To show the flexibility of the particle 
algorithm with respect to a choice of a time integration scheme, we chose a symplectic 
integrator of fourth order accuracy (the PEFRL algorithm of Omelyan et al. [46]). For a 
choice of time step At = 0.01, energy conservation is virtually perfect, at approximately 
j^Q-is ijiaximum relative error. 
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Figure 4: Linear growth and saturation of t he fi rst fo ur h armonics in the beam— plasma problem com- 
puted using potential-based particle model, | |32| l and | |33| l, for ni,/no = 10"'' and the same simulation 
parameters as in Fig. [2] 



Similar results, shown in Fig. [4j are obtained when the beam-plasma instability sim- 
ulation is performed with the potential-based model, (32) and (33). (The regions where 
the growth rates are determined are indicated by dots and correspond to the same regions 
used in Fig. [2j) The equations of motion were integrated with a second-order Runge- 
Kutta method with time step At = 0.001. No time-splitting was used, i.e., all particle 
and field data are known at common points in time. The number of grid points was 2048, 
the number of particles per cell was 4, pk was cubic in ^a, corresponding to the shape 
Si (see Table A.l ). The growth rate of the first few harmonics is in excellent agreement 
with the truncated Fourier series model. Energy conservation for this model is also very 
good, with relative energy error of less that 0.6% (not shown). 

The examples of Figs. [2]-[4] demonstrate that energy conserving algorithms perform 
reliably in this benchmarking test and have low noise due to the freedom to choose 
smooth particle shapes. 
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(d) 5, (Quartic) 
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Figure 5: Dependence of the fc = 1 growth rate in the beam— plasma problem on grid size for the potential- 
based particle model, l |32| and ( |33[ l. Plotted is the difference between the growth rate calculated using 
a given value of Ax compared to the growth rate computed with very high spatial resolution for various 
values of At. As the plot shows, the method is second-order, as expected by our approximation of C. 
The panels are labeled with the particle shape S and the resulting order of (see Table|A.l|. 



Figure [5] shows the dependence of the fc = 1 growth rate of the beam-plasma problem 
on the spatial resolution in the potential-based model (|32| and (33). Plotted is the 



difference between the growth rate calculated using a given value of Aa; compared to 
the growth rate computed with very high spatial resolution (and the specified value of 
At.) As expected from our approximation to C, growth rate is second order in the grid 
spacing, Ax, regardless of the particle shape or time-step. 

The next example illustrates an important property of these energy conserving algo- 
rithms, namely, that energy conservation depends solely on the properties of the time 
discretization. We consider a linear plasma oscillation, with electric field amplitude of 0.1 
and integrate the equations of motion with a second-order Runge-Kutta method. The 
relative energy error at t = 400 is plotted against the time step At for various particle 
shapes (see Table A.ll and grid resolutions. Figures [6] and [7 show the relative error 
for the potential-based, (32) and (33), and field-based, ( |104 1-( 106 1, energy conserving 
algorithms, respectively. The scaling with time step for all particle shapes is ~ 0{At^); 
the exception is only for the potential-based method with linear pi^. In this case the 
force on a particle is discontinuous, i.e., it has jumps as a particle moves from one cell to 
another. Therefore, for the potential-based formulation linear particles are not recom- 
mended. Interestingly, in the case of linear particles in the field-based formulation, this 
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Figure 6: Conservation of energy for tlie potential-based particle model, (j32| and ( |33[ l. The relative 
energy error is shown as a function of At for various spatial resolutions and particle shapes. As expected, 
the energy error depends only on the temporal discretization. The equations of motion were solved with a 
second order Runge-Kutta integrator. The panels are labeled with the particle shape S and the resulting 
order of (see Table [ATT| | . 
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(a) Sx (Linear) 
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Figure 7: Conservation of energy for the field-based particle model, | |104| l-( [l06| l). The relative energy 
error is shown as a function of At for various spatial resolutions and particle shapes. As expected, the 
energy error depends only on the temporal discretization. The equations of motion were solved with a 
second order Runge-Kutta integrator. The panels are labeled with the particle shape S and the resulting 
order of (see Table [ATT| | . 
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Figure 8: Conservation of energy for a standard (momentum conserving) PIC algorithm. The relative 
energy error is shown as a function of At for various spatial resolutions and particle shapes. The 
equations of motion were solved with a second order Runge-Kutta integrator. The panels are labeled 
with the order of the charge-deposition/force-interpolation spline used; see Ref.[21j. 
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deficiency does not sliow as strongly and tlie scaling has the same trend as for smoother 
particle shapes, Fig. [7] Simulations with the truncated Fourier basis particle model ex- 
hibit similar behavior (not shown). Figure [s] shows energy conservation in the standard 
(momentum conserving) PIC algorithm. As expected, the relative energy error in these 
algorithms depends on both the time step as well as the spatial grid resolution. Note that 
for the same particle smoothness and time step, the energy conserving algorithms have a 
(much) smaller relative energy error for At < 10~^; for smaller At, energy conservation 
in the PIC algorithm is limited by the maximum number of grid points being fixed at 
512. These examples combined with the example in Fig. |5] demonstrate that our method 
has overall accuracy of second order in both time and space. 

5. Conclusions 

We have derived time-explicit, energy-conserving algorithms based on two approaches: 
Lagrangian in terms of potentials and a Hamiltonian with a non-canonical Poisson bracket 
in terms of fields. These models are derived without specifying any particular spatial 
or time discretization scheme, accuracy, or particle shape. Our general method allows 
the Lagrangian-based derivation to relax a number of restrictions imposed previously. 
Continuous quantities are reduced by performing either a grid reduction (i.e., finite dif- 
ferences) or truncated bases. When a grid reduction is used, mass matrices do not 
appear, which decreases computational load and improves the efficiency of memory us- 
age. The important role of the particle shape and its relation to force interpolation is 
exhibited. A relaxed choice of particle shape helps decrease numerical noise in energy- 
conserving algorithms. A Hamiltonian derivation is presented here for the first time. 
The method uses a reduction of both the Hamiltonian and the non-canonical Poisson 
bracket. Since its formulation is in terms of fields, it avoids solving Poisson's equation. 
A model conserving both energy and momentum is derived and the conditions that make 
this possible are described. Its derivation uses the relation between conservation laws 
and Lagrangian symmetries, thus emphasizing the power of variational principles. Nu- 
merical benchmarking confirms the improvements in our algorithms. It is shown that 
conservation of energy in all particle models derived here only depends on the accuracy 
of time integration; in comparison, energy conservation in PIC depends on both the grid 
spacing and time integration accuracy. We restricted our discussion to the case of a 
one-dimensional, nonrelativistic, unmagnetized, electrostatic plasma. The generalization 
to three dimensional, relativistic, electromagnetic plasmas is straightforward and will be 
presented elsewhere. It is shown how to increase overall accuracy (in space and time) 
beyond second order. 
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Appendix A. Particle Shapes 



We present charge-deposition rules Pi{(,a) for a variety of particle shapes based on 



(23) where the interpolation uses linear finite-elements; see Table A.l and Figure A. 9 



The shapes S'o, and 5*2 correspond to the usual PIC particle shapes scaled by l/h. 
(The standard PIC definition normalizes the particle shape to have area h whereas we 
normalize our shapes to unity.) This helps to explain the poor energy conservation 
observed with linear deposition (see Figures |6] and [s]) ; such deposition correspond delta- 
function macro-particles. As a result, any deposition scheme used should be at least 
quadratic in (*-e., at least C^). All of these charge deposition rules are second- 
order accurate; higher-order accuracy is obtained by with a correspondingly higher order 
interpolation method, e.g. using quadratic finite-elements would lead to a third order 
accurate interpolation. 




Figure A. 9: Particle shapes and corresponding charge deposition. Charge deposition corresponding to 
various particle shapes S{x — ^ci). For a given particle position ^a, is the nearest grid point. We show 
the entire range of for which is non-zero so that the effective particle shape is apparent. The range 
of for which — 5a| < h/2 is indicated by the dashed grey lines. See Table [aTT] for the definitions 
of 5„. 

As discussed in Section [2] in our formulation, there is no requirement for the particle 
size to be tied to the grid spacing. To demonstrate this freedom, consider 

2 

S{x^U = j;{ (A.l) 



. otherwise. 
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From Eq. ( 23 ) , using linear finite elements for interpolation, and with the same definition 
of A as above [A = (^q, — Xk)/h], we have for A > 



and for A < 



Pk 



Pk+l 



Pk = 



Pk+l — 



Pk-1 = 



5 

6 
1 
12 

1 
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-2A2 
1 . 



2 A' 



A 



A^ 



3 ' 



A' 



A' 



A' 



(A.2) 



(A.3) 



This particle shape results in a cubic deposition scheme (i.e. the pk are in 



involving only three grid points (or 2 cells); see Figure A. 10 In contrast, the usual PIC 



cubic deposition involves four grid-points corresponding to a particle three cells in extent. 
Because of the linear finite elements used, this still produces second order accurate force 
interpolation. 



hS{x-U 




Pk iL ) 




Figure A. 10: Charge deposition corresponding to various particle shapes S{x — ^a)- For a given particle 
position ^a, is the nearest grid point. We show the entire range of for which pj. is non-zero so 
that the effective particle shape is apparent. The range of for which — < h/2 is indicated by 
the dashed grey lines. See ( |A.10[ | for the definitions of Sn- 



Appendix B. A fluid kinetic hybrid model 

Under certain conditions, plasma electrons, plasma ions, or both, may be well ap- 
proximated by (charged) fluids. Typically in hybrid models, one species is described as 
a fluid while the other as particles (i.e., kinetically). Here we describe a single species 
with a fluid-kinetic hybrid and make no assumptions about the inter-mixing of the fluid 

32 
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Pk+2 = ^ (A + i) 
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96 ' 24 



|A^ + iA^ 



+ iA+lA^ 



1A3 



^A* 



Table A.l: Charge deposition corresponding to various particle shapes S. For a given particle position 
^a, let k be the nearest grid point and A = {^cx — Xk)/h. All pi other than those listed are identically 
zero. 
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and kinetic elements. A prototypical system is a low-charge electron beam propagating 
in a cold, quasi-neutral plasma. Only the beam {i.e. the tail of the distribution) needs 
to be treated in a fully kinetic manner. The bulk plasma can be represented as a fluid. 
When the bulk plasma thermal velocity is small compared to both the velocity of the 
electron beam and the phase velocities of plasma waves arising from the beam-plasma 
interaction, we may take the fluid to be cold. 

It is computationally advantageous to use such splitting when the kinetic population 
of the plasma is small compared to its fluid-like population. The numerical noise in 
the hybrid description can be much lower compared to the kinetic description of the 
entire plasma. In addition, there can be a large computational speedup due to using a 
fluid description for the larger fraction of the plasma since only time advance of gridded 
quantities is required. 

For concreteness, we assume stationary ions and mobile electrons. The cold electron 
fluid distribution function is approximated as a delta function in velocity space 

/(x, V, t) = 7i(x, t) 5[v - A x(0] , (B.l) 

where n(x, t) is the fluid density, DtK{t) is the velocity of a Lagrangian fluid element 
{Dt = dt + Uid/dxi being the convectional derivative), and u is the Eulerian fluid 
velocity. We express the fluid velocity in terms of the velocity potential T and the 
Clebsch variables a and /3 as u = VT -|- aV/3 with V x u = Va x V/?. The kinetic 
electron distribution is given by expression ([t]) and the complete distribution function is 

/(x, V, i) = n(x, t)8[w - Dt x(t)] + ^ u;„5[x - ^„ [t)] ^ [v - {t)] . (B.2) 

a=l 



The reduction proceeds as before, by substitution of (B.2 1 into Low's Lagrangian, ([3|. 



Since the Lagrangian is linear with respect to the distribution function we can con- 
sider separately the fluid and kinetic contributions to / simply adding the resulting 
Lagrangians. The reduction of the kinetic contribution is identical to Section [2] and the 
reduced Lagrangian (without the field contribution) is given by 

Np Np Ng 

^^Wail- qs^^WaPi{^a)'Pz- (B-3) 
a— 1 a— 1 i—1 

Below we consider only the fluid and field contributions ^ together. At the end we add 
the particle, fluid, and fleld contributions. 

Again for simplicity, we specialize to a non-relativistic, electrostatic system. One 
may formulate the fluid variational principle in Eulerian picture using the fluid density, 
velocity potential, and the Clebsch variables as independent variables [17] : 



£f = —TLs / d xn 



^(VT + aV/3)' + t + a/^ 



1 

8^ 



d^x{W(pf-qs Id-^xnip. (B.4) 



The second (field) term in (B.4) must be included only once in the final Lagrangian of 
particles and fiuid. As a further simplification, we develop the fiuid model in one spatial 
dimension; we also assume that ions form a uniform (immobile) background with density 
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no- In one spatial dimension, a = (3 = and the fluid experiences electrostatic force 
through the electrostatic potential tp. The Lagrangian for this case is 



dx 'I 



+ ^ jdx {S/f f-qs jdx nip-q^n^'"''^ j dx ip . (B.5) 



To proceed with the reduction, we first choose a linear finite element basis and expand 
all dependent variables. We will see, this approach leads to mass matrices, which need 
inversion. Then we use grid-based reduction, which eliminates mass matrices altogether, 
and therefore such numerical model has computational advantage. (Other choices, such 
as a truncated Fourier basis, may be more appropriate for some applications; nonetheless, 
the derivation proceeds along similar lines.) 
The fluid variables are represented as 

Wo 



T(x,t)=^T,(t)vI/<^>(x), 

i=l 
i=l 



(B.6) 



Substitution of (B.6) into the fluid Lagrangian (B.5 1 gives 



1 

m, 

2 



'ijifiipj - gs ^ Mij ipi n.j - Qih n''™' ^ (pi, (B.7) 



where h is the grid spacing. 



Fijk = Idx 



dx 



dx 



(B. 



and Kij , and Mij are given by ( 26 ) and ( 99 1 , respectively. (We note that the coefficients 
Fijk are symmetric in i and j but not symmetric in all three indices.) 

This is a finite degree-of-freedom Lagrangian and thus requiring the action to be 
stationary, leads to the usual Euler-Lagrange equations 



-0, 

= 0, (B.9) 
0. 



d dCp 




dtdTk 
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dt dipk 
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dCp 


dt dfik 
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The Euler-Lagrange equation for ni yields the cold fluid momentum equation: 

N„ 



1 



FijiTiTj - ^ Miiti - Mil Vi = 



or 



JL Z = „ 

m, 2 



The Euler-Lagrange equation for T/ gives the fluid continuity equation: 



Ng Ng 

rus E Mijhj + TOs ^ FijkTjUk = 



i=l 



or 



(B.fO) 



(B.ll) 



(B.12) 



(B.13) 



Finally, the Euler-Lagrange equation for ipk gives Poisson's equation (without the particle 
contribution): 



(B.14) 



Eqs. (B.11)-(B.14) are a complete set of equation for a fluid description of a plasma. 
To obtain the hybrid particle-fluid model, these fluid equations must be supplemented 
by the particle equation of motion and the particle contribution to Poisson's equation: 



Co 



1^^. 



(B.15) 



Na 



i=l 

^K.j^j = -47r^^u;,p,(e«)-47r^^Af,,nj-4^g:n"™', (B.16) 



a=l 



where Pi{S,a) is given by (23). The complete set of hybrid fluid-particle equations is given 
by Eqs. ( |BJ1] ), ( |rI3| , ( |R15| ), and (|b7T6|. The conserved energy is 



Wr. 



- ^ Wai'i ~ 2 ™^ ^ -PjjfcTi Tj nk - Kijipiipj. (B.17) 



Q = l 



i,j,k=l 



We see the appearance of the inverse of the mass matrix in Eqs. (B.ll) and (B.13). 
Having to keep the (dense) inverse of the mass matrix may consume too much computer 
memory and make computation times longer. Therefore, we show that with the use 
of the grid-based reduction mass matrices to not appear and no matrix inversion is 
necessary. The grid-based reduction simply tells us to use a numerical integration rule 
and finite differences to reduce the integrals of continuous quantities. Thus, all bi-linear 
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combinations of continuous quantities reduce to sums over the grid index. With these 



rules, the Lagrangian (B.7) becomes 



with 



ij",fe=l 



N„ 



Stt 



ijipiipj ~ Qsh^^ipi Tii - qihn'""''' ^ ipi, (B.18) 



ijk 



hDijDik- 



(B.19) 



Kij was defined in ( 29 ) and the finite differencing operator Dij may be chosen as the 
second order accurate centered differencing Dij = ((^^.i+i — Si i_i)/2h. It is now clear 
how to obtain the corresponding modif ied eq uations of motion: replace Mij by h (and 
M^by /i-i) and use the coefficients te.l9| in place of te.8|. The energy expression 



(B.I7I must also be modified accordingly, ibtal momentum for this hybrid model is not 



conserved due to the use of grid. 
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